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Self Organized Criticality in a two dimensional Cellular Automaton model of a magnetic 

flux tube with background flow 
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Abstract 

We investigate the transition to Self Organized Criticality in a two-dimensional model of a flux tube 
with a background flow. The magnetic induction equation, represented by a partial differential equation 
with a stochastic source term, is discretized and implemented on a two dimensional cellular automaton. 

The energy released by the automaton during one relaxation event is the magnetic energy. As a result 
of the simulations we obtain the time evolution of the energy release, of the system control parameter, 
of the event lifetime distribution and of the event size distribution, respectively, and we establish that 
a Self Organized Critical state is indeed reached by the system. Moreover, energetic initial impulses in 
the magnetohydrodynamic flow can lead to one dimensional signatures in the magnetic two dimensional 
system, once the Self Organized Critical regime is established. The applications of the model for the 
study of Gamma Ray Bursts is briefly considered, and it is shown that some astrophysical parameters 
of the bursts, like the light curves, the maximum released energy, and the number of peaks in the light 
curve can be reproduced and explained, at least on a qualitative level, by working in a framework in 
which the systems settles in a Self Organized Critical state via magnetic reconnection processes in the 
magnetized Gamma Ray Burst fireball. 

Keywords: Magnetic reconnection; Self-Organized Criticality; MHD (magnetohydrodynamics); Stars: 

Gamma-ray bursts. 


1 Introduction 

Many important natural dynamical systems show the presence of long-range spatial and temporal correla¬ 
tions. For example, spatial scale-invariance is observed in fractal geographical and topographical structures - 
mountain ranges, river basins, etc. [28) . while temporal scale-invariance in the form of l//-like power spectra 
is observed in such diverse phenomena as star flickers and earthquakes [351 f] [2] [TO] . The spatial scale in¬ 
variance in particular is in marked contrast to the typical behaviour of equilibrium thermodynamic systems, 
where such an invariance can be realized by tuning a parameter (e.g., the temperature) to a critical value. 
To describe scale-invariance that occurs in dynamical systems without the explicit timing of a parameter, a 
model pioneered by [23] and later popularized as Self-Organized Criticality (SOC) by [3] and [4], has begun 
to develop. 

As a paradigm of self-organized critical behaviour, [3j introduced the sandpile model, which they im¬ 
plement with a Cellular Automaton (CA). A one-dimensional simple local version of the CA can be briefly 
described as follows 12315]. We start by prescribing two positive integers, z c and n. Then to each site 
i = 1, L, we assign an integer h(i) > 0, representing the height of sand at location i. The slope of the 
sandpile at i is given by z(i) = h(i) — h(i + 1). The system is assumed to be closed at its left boundary 
( i = 0), and open at the right one (i = L + 1). The dynamical process consists in dropping a single grain of 
sand onto a randomly chosen site. If the slope at any site i exceeds z c , then n grains of sand fall from i onto 
the site i + 1. This process may cause the slope to exceed 2 C at adjacent sites. If necessary, n grains fall to 
the right from any site at which the slope now exceeds z c . The process continues until all slopes are less than 
z c . This may be accomplished either with or without the loss of sand from the right boundary. The entire 
event caused by dropping a single grain is called an avalanche. The number of sites from which sand has 
fallen is called the size of the avalanche. After an avalanche, a grain is dropped again onto a randomly chosen 
site, hence another avalanche occurs, and the process continues indefinitely. The obtained numerical results 
nna indicate a power law distribution in the sizes of avalanches. This is exactly the type of distribution 
found at the critical point of traditional equilibrium systems. Quantities related to the avalanche sizes also 
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have power law behaviour. Thereby one can define a set of critical exponents and scaling relations for the 
system [5j. 

As an alternative to the numerical study of SOC in cellular automata, one can consider models based 
on differential equations that describe the fluctuations of a conserved quantity m ns ii mi CEB- More 
specifically, systems described by diffusion equations driven by noise terms shows characteristics similar to 
SOC, with the numerical solutions of partial differential equations on finite grids yielding cellular automata 
with real-valued states. For example, the basic one dimensional diffusion equation du/dt = (l/2)d 2 u/dx 2 
can be transformed into the discrete equation u(x,t + 1) = (1/2) [u(x + l,i) + u(x — l,t)] [35] ■ Conversely, 
the diffusion limits of suitable cellular automaton models yield partial differential equations. The solutions of 
such diffusion equations exhibit scale-invariant behaviour, similar to that observed in sandpiles and related 
discrete models. It turns out that SOC can already be understood on the level of the linear stochastic 
Langevin equation jS], This is due to the fact that algebraic decay of spatial correlations is a reflection in 
spatial dimensions of the algebraic decay in time for systems with conservative dynamics. 

The possibility that SOC results in some specific astrophysical phenomena has been intensively inves¬ 
tigated recently, and the study of SOC has become a major field of research in astrophysics mum. In 
particular, an interesting model to describe the properties of the accretion disks around black holes was 
introduced in :3II;, where it was suggested that the inner portions of black hole accretion disks may be 
in a self-organized critical state. Consequently, l//-like X-ray fluctuations are produced, in spite of random 
mass input. Viscous diffusion processes are also incorporated in the model. Hence in this approach mass 
accretion occurs either by an avalanche, which is triggered when the mass density of the disk exceeds some 
critical value, or by gradual gas diffusion. This initial model was further modified in [36] . where a gradual 
diffusion which occurs regardless of the critical condition was introduced. In m it was shown that with 
only one free model parameter, the cellular automaton model can reproduce some observational data of Cyg 
X-l. The relativistic effects were included in the model in [41] . where it was shown that the CA model can 
produce light-curves and power-spectra for the variability that agree with the range observed in optical and 
X-ray studies of AGN and X-ray binaries. It was also pointed out that when general relativistic effects are 
incorporated, important differences do appear if the disk is viewed from directions far from the accretion 
disk axis. The magnetic activity of an accretion disc using a probabilistic cellular automaton model was sim¬ 
ulated in [32] . The model is based on three free parameters, the probabilities of spontaneous and stimulated 
generation of magnetic flux above the surface of the disc, and the probability of diffusive disappearance of 
flux below the surface. This approach allows steady accretion in a disc by the action of coronal magnetic flux 
tubes alone, and if the effective viscosity caused by coronal loops is expressed in the usual Shakura-Sunyaev 
alpha parameter of viscosity, one can obtain numerical values which are in good agreement with observations. 
A modified SOC model based on cellular automaton mechanism for producing lognormal flux distribution 
was presented in [24]. In this model the energy released in the avalanche and diffusion in the accretion disk is 
not entirely emitted instantaneously, with some part of the energy kept in the disk. Thus the disk increases 
its energy content, so that the next avalanche will be in higher energy condition, and more energy will be 
released. Hence the later an avalanche occurs, the more amount of energy is emitted from the disk. 

On the other hand it was found observationally that a class of Gamma-Ray Bursts (GRBs) called Soft 
Gamma Ray Repeaters (SGR) do exhibit SOC characteristics [21 OS] • From the statistics of the SGR 1806- 
20 bursts it was shown that the fluence distribution of bursts observed with different instruments is well 
described by power laws with indices 1.43, 1.76 and 1.67, respectively. Thus, it turns out that the SGR 1806- 
20 bursts behave in a self-organized critical manner, similarly to earthquakes and solar flares m- From an 
astrophysical point of view this behavior suggests that the energy sources for SGR bursts are crust-quakes, 
due to the evolving, strong magnetic field of the neutron star, rather than any accretion or nuclear reaction 
processes. A very similar result was obtained for the case of the X-ray flares of GRBs with known redshifts by 
[39] , who have shown that X-ray flares and solar flares share in common three statistical properties: power- 
law frequency distributions for energies, durations, and waiting times. All these distributions are specific 
for the physical framework of a SOC system [8J. Hence we have the interesting result that the statistical 
properties of X-ray flares of GRBs are similar to solar flares, and therefore both can be attributed to a 
SOC process. As suggested in [39], both types of flares may be driven by a one dimensional SOC magnetic 
reconnection process. On the other hand the X-ray flares of GRBs cold be produced astrophysically in 
ultra-strongly magnetized millisecond pulsars, or long-term hyperaccreting disks around stellar-mass black 
holes. 
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The possibility of the one dimensional SOC modelling of X-ray GRB afterglows, as well as the possibility 
of appearance of Self-Organized Criticality in an one dimensional magnetized flow was carefully investigated 
i n [151 - A simplified one dimensional grid was used to model the evolution of the magnetized plasma 
flow. Diffusion laws similar to those used to model magnetic reconnection with Cellular Automata in various 
astrophysical phenomena were implemented in the model, as well as a background flow. Under the assumption 
that the parameter relevant for X-ray afterglows is the magnetic field, the magnetic energy released by one 
volume during one individual relaxation event was computed. The obtained results show that indeed in this 
system SOC is established. The possible applications of this model to the study of the X-ray afterglows of 
GRBs was also briefly considered. 

However, from a strictly theoretical point of view one dimensional magnetic reconnection is not possible. 
This raises the interesting and important question of clarifying the relation between GRB observations, 
showing that indeed one dimensional SOC occurs in the astrophysical process, and their possible relation 
with one dimensional magnetic reconnection models. In order to answer to this question in the present paper 
we generalize the one dimensional model introduced in [ T8] . 

Hence in order to describe the temporal evolution of GRBs we develop a two-dimensional CA simula¬ 
tion algorithm, which implements the magnetic induction equation in the MagnetoHydrodynamic (MHD) 
approximation framework. The bulk advection motion is, as a novelty with respect to other models, taken 
explicitly into account for various velocity profiles. The end purpose is to find a model which is simple to 
implement but which still captures the important macroscopic characteristic of observed X-Ray afterglows 
in GRBs. The main result of our analysis is that even the rigorously implemented physics is 2D, energetic 
initial impulses in the MHD flow lead to ID signatures once the SOC regime is established. 

The present paper is organized as follows. The simulation procedure is described in Section [2j and the 
discretized MHD equation for the magnetic induction are written down. The simulation results, and some 
of their astrophysical implications, are presented in Section [3] We discuss and conclude our findings in 
Section [4j 

2 Setup of the simulation 

Although our goal is to keep the presentation in this Section self-consistent, the reader is referred to [18] for 
technical details of the simulation that might be absent from our current presentation. In the following we 
concentrate on the evolution of the magnetic field in the cosmic environment. From a physical point of view, 
under certain conditions, the magnetic field topology can change suddenly, thus exhibiting a transition in 
its behaviour. The transition can be understood and quantified in terms of the Laplacian of the magnetic 
field achieving a critical value |8j. The evolution equation for the magnetic field is the magnetic induction 
equation. It describes the time and space evolution of a magnetic field in a non-ideal plasma medium. The 
magnetic induction equation is derived under the assumption that the MHD approximation holds in the 
plasma [33], and that there is a background flow. We will present below a cellular automata approach to 
plasma dynamics. The information about the relevant physical parameter describing the dynamics of the 
system is stored in the cells of the grid [7]. 

2.1 The magnetic induction equation 

We begin our presentation of the mathematical model to be simulated with the induction equation j3D] [33] 

— = V x (Fx BJ + ? ? V 2 H, (1) 

where B = B(x,y,z,t) is the magnetic field, v = v(x,y,z,t) is the plasma velocity, and p is the total 
magnetic diffusivity coefficient, ?/ = l//xocr, where a is the electrical conductivity, and po is the magnetic 
permeability in a vacuum. We define the control parameter in the plasma flow as 

G = ~V 2 H, (2) 
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to be computed for the four neighbours in a two dimensional grid (left, right, up, down). The configuration 
we are working with is 

B = (B x (y,z,t), 0,0), (3) 

v={0,0 ,v z (z,t)). (4) 

The magnetic field evolution equations then become 


dB x d{v z B x ) | d 2 B x t d 2 B x 

~df = + + 


G 


X 


1 fd 2 B x d 2 B x \ 

4 ^ dz 2 + dy 2 ) ' 


( 6 ) 


In the numerical simulations, the time dynamics of the system is usually implemented as one of two 
different regimes (advection or diffusion), as a function of the value of the control parameter. A connection 
to the physical process behind this switch between one regime or another can be made by the following 
reasoning. In astroplrysical conditions the classical resistivity is very small, and the magnetic field behaves 
macroscopically as if the diffusion term in the magnetic induction equation would be zero. This behaviour 
is controlled by the magnetic Reynolds number R m , 
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UL 
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(7) 


where U and L are a characteristic velocity and a space-scale, respectively. From a physical point of view 
R m is a measure of the size of the advection term, V x (y x B ^, as compared with the size of the diffusion 

term, ?;A B. When R m >> 1, the diffusion term in Eq. (JT|) is negligible, while for R m « 1, the evolution 
of the magnetic field is purely diffusive. 

However, under certain conditions, in relatively small volumes, the diffusive behaviour becomes dominant, 
and the magnetic field lines reconnect. To model this process in our grid, we assume that the astropliysical 
evolution occurs in two different regimes, with the switch between these two regimes determined by the 
behaviour of the control parameter. The control parameter tells us what is the value of the difference 
between the magnetic field between one grid point, and its neighbours. By definition, this control parameter 
is thus local, and its characteristic scale l is very small. If L is the characteristic length scale of the simulation, 
then in our case l/L is at least 0.002. If the control parameter exceeds a certain threshold, then what happens 
locally becomes worthwhile inspecting. Since we may assume that the velocity does not change in order of 
magnitude, and since the diffusivity y is constant (the plasma properties do not change), the ratio between 
the macroscopic Reynolds number and that of the local Reynolds number is of the same order of l/L. This 
can be viewed as a reason why the control parameter changes the behaviour of the Reynolds number. 

To summarise, the advective regime is the main framework in which we develop our model. If the control 
parameter becomes critical the magnetic field evolution is given, for a brief period of time, by a diffusive 
behaviour. Once this local criticality is relaxed, the control is given back to the advective evolution. 

Analytically, 

dB x _ j - d{B d X z Vz) . hi S h R m » 1, 

dt j ?? + 7 ^) , low R rn « 1. 

We take the term as the stochastic source term and denote it by S(y, z, t). 

The evolution equations are brought to dimensionless form by the scalings 


t = aT, B x = b(T, Z , Y)B 0 , z = /3Z, y = yY, 

v = v 0 V,G x = gB 0 //3 2 , (9) 


where a, Bq, /3, 7 and vq are constants. With these transformations, the system description becomes 

Q2 b _ fl2 


db d 2 b 
fff^ k 'dZ I + ak df' 


or? ^ 
k 7 2 > 


( 10 ) 
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for the diffusive behaviour, and 
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( 11 ) 


for the advective behaviour, where s(Y,Z,T) = aS(y 7 z,t)/B 0 ; the dimensionless expression for the control 
parameter is 

1 ( d 2 b d 2 b\ 

9 ~ ~4 \dZ 2 +a dY 2 ) ' ( ^ 

In the Solar Corona the magnetic diffusivity 77 can be approximated as rj « 10 4 (T/10 6 K) 3 ^ 2 cm 2 /s, 
while the diffusion region has a size of L ss 2 x 10 s cm [33]. By adopting a characteristic velocity of the order 
of U = 10 ' 3 cm/s, it turns out that the magnetic Reynolds number of a plasma with temperature T = 10 6 K 
is around R rn = 2x10' >> 1. In this case the diffusion term in the magnetic diffusion equation Eq. Ill is 
negligible. However, in Solar Physics there are several important exceptions, like, for example, in the neigh¬ 
bourhood of magnetic neutral points and lines, during magnetic reconnection, and during solar flares, when 
the diffusion processes play an important role in the understanding of the corresponding processes |34j . In the 
case of the Gamma Ray Bursts, one can estimate the magnetic Reynolds number by adopting for the magnetic 
diffusion coefficient rj the perpendicular resistivity in a strong magnetic field r]± = 1.3 x 10 13 x Z In A/T 3 / 2 
cm 2 /s, where A = 3/2e 3 ^/k 3 T 3 /irn is the Coulomb logarithm, with k denoting Boltzmann’s constant and 
n the particle number density, respectively [42| . By assuming physical conditions specific for GRBs, we 
find r)± ss 1 cm 2 /s, leading to a large magnetic Reynolds number of the order of R m ~ 10 24 . However, by 
assuming for the magnetic diffusivity in GRB magnetized plasma its maximum value r/B, corresponding to 
the Bohrn diffusion, we obtain r]B ~ r^c = 7 e m e c 3 /eB 7 where re = 7 e m e c 2 /eB is the co-moving frame 
cyclotron radius, and is the electron Lorentz factor. Hence for the magnetic Reynolds number in the 
Bohrn diffusion approximation f? m ,Bohm we obtain R m ,Bohm ~ 3.4 x 10 12 7 f r 1 L^ // 2 2 r^ 2 -\/(T/(l + a) cm 2 /s, 
where L Wt 52 is the luminosity of the GRB ejecta in units of 10 52 erg/s, T 2.5 is Lorentz factor of the ejecta in 
units of 10 2 ' 5 , and a is the magnetization parameter [42] . Hence for high values of the particle and ejected 
wind Lorentz factors, and for relatively low Gamma Ray Burst luminosities, in the Bohm diffusion limit 
considerably small magnetic Reynolds numbers may characterize the plasma ejected by the GRB explosion. 
For example, for y e = 10 10 , L Wt 52 = 1, T 2.5 = 100, and a » 1, -R m ,Bohm ~ 3.4 x 10 -2 . Therefore, during 
the post-explosion expansion of the relativistic fireball in a strong magnetic field both very high and low 
magnetic Reynolds number regimes may be present, and this possibility must be taken into account when 
analysing the dynamical behaviour of the GRB lightcurves. 


2.2 The discretized model 

Following the regular discretisation procedure (with indices i, j, k standing for the Y, Z and T coordinates), 
we obtain for the above equations 
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(A zy 


\bi,j+i,k + bi,j—i,k 2bij t k] 


(13) 
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for the critical parameter. 

We will consider flows with a velocity decreasing as the spatial grid index increases, such that 




( 16 ) 


and, even more, propagating flows such that the velocity is zero for points not yet reached by the wave-front 
along the vertical direction Z, i.e., 


For the behaviour of the wave-front with respect to the Y axis we consider a Gaussian function centred on 
the current value of the Z coordinate, of a width fixed in the code, but such that the integral of the Gaussian 
with respect to Y is 1. 

When the critical threshold has been reached in a point, the diffusive behaviour of the induction equation 
is 

4 

bi,j,k -FI t bi.j,k (18) 


and the redistribution is 


bi,jz Ll,fc+1 t bi jk -f 

a 


(19) 

( 20 ) 


The energy released by each volume during an individual relaxation event is the magnetic energy lost in that 
volume 

e r= f [I [ b2 x (in) - B 2 X (out) ] dxdydz (21) 

Z ^0 Jdz Jdy Jdx 


and in dimensionless form it is 


£r = 


= 7T~ f [ [B 2 * (in) - B 2 X (out)} dYdZ , 

2/^0 JdZ JdY 


( 22 ) 


where 


Hence we find 


Er = eu 


hBl 
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€-k 


EE 


8 _ 16 2 

^Oij t k*gi,j,k* 25 '^ i ’ 



(23) 

(24) 


The star on the timestep counter k represents the fact that within one time step k of the simulation, the 
same cell, due to next neighbour interaction might become unstable more than once. The k* is a subdivision 
of the simulation time step and it is non-zero as long as the cell is unstable. 


2.3 The simulation procedure 

The simulation procedure is described below as follows: 

1. Initialization: a two dimensional grid with Ny x Nz cells is initialized in each cell with the value bo ; 
the initial flow velocity (Vo,o,o) is some multiple of the characteristic Alfven speed for the configuration, 
given by y; 

2. Evolution: for each k th timestep in the interval 1, Nt, the evolution of the system is as follows: 

• Since the upward flow with velocity V is deterministic, one can formally know what cell jh the 
flow has reached at the time step k. For fixed current Z coordinate, the Y cells associated to it are 
updated with values drawn from a Gaussian centred on the current Z , with standard deviation 
Ny/ 20 (advective behaviour); 
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Figure 1: Left panel: plot of the Counts vs. time data of GRB 140919636 on 2014-09-19, obtained with 
the FERMI Gamma-ray Space Telescope with detector 0 dzj. Right panel: Simulated energy release as a 
function of time, for x = 50 and the 14 =const. velocity profile, caused by a single GRB pulse impinging at 
the base of the simulation grid. 


• A random pair of numbers {ki, kj}, ki £ 1 ,Ny, kj € 1 ,Nz is chosen and updated as 

Si,j,k = h,j,k{ 1 + A Te), (25) 

where e < 1 is a positive small number (stochastic loading); 

• The mean value of the magnetic field b is calculated; each component of the grid is then scaled 
with respect to this mean value. The critical parameter g cr is taken as 10% of this value. 

• Scanning: The control parameter is calculated with Equation (1161) for each cell in the grid and 
stored into a matrix. If any one cell has an absolute value higher than the control parameter g cr 
a flag is triggered; 

• Redistribution: If the flag has been triggered, the control parameter matrix will be searched 
for positive values that are higher than g cr . If any are found, part of the cell’s content will be 
redistributed to its nearest neighbours according to Eqs. m-m- If no positive values were 
found, the matrix will be searched for negative values with absolute value higher than g cr . This 
sweep is done while cells with value higher than the critical parameter are found. A variable A4 
stores the number of such events for each step k. Another variable Sk stores the number of cells 
reached by the critical flow. 

3. Results: A4 represents the number of events needed to fully relax the grid at each time step k; the 
vector N is used to produce the lifetime distribution, D(N ); similarly, Sk represents the number of 
cells reached by the critical flow at each time step k: the vector Sk is used to produce the lifetime 
distribution, p(S). 

The connection between grid parameters, observational parameters and simulation output is given in 
Table [T| 

In order to obtain a description of the stochastic behaviour of the light curves we use the slope of the Power 
Spectral Distribution (PSD) of the luminosity. The slopes of the PSD curves can provide some important 
insights into the nature of the physical mechanisms leading to the observed variability of the astropliysical 
source. If A is a fluctuating stationary physical quantity, with mean fix and variance cry, respectively, then 
we define the autocorrelation function for A' as follows [551 

o / ((-^« — A 1 A') (A s+r — /i)) ,- 0 c\ 

Rx{t) = -2-> ( 26 ) 

a x 
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Figure 2: Left panel: plot of the Counts vs. time data of GRB 080805496 on 2008-08-05 obtained with 
the FERMI Gamma-ray Space Telescope with detector 0 [15]. Right panel: Simulated energy release as a 
function of time, for x = 10 and the 14 =const. velocity profile, caused by a single GRB pulse impinging at 
the base of the simulation grid. 




Figure 3: Left panel: plot of the Counts vs time data of GRB 101126198 on 2010-11-26, obtained with the 
FERMI Gamma-ray Space Telescope with detector 7 ild . Right panel: Energy release as a function of 
time, for \ = 10 and the l/Vk velocity profile, caused by a single GRB pulse impinging at the base of the 
simulation grid. 


Table 1: Parameter correlations. The Dimensionless column contains parameters which are set beforehand, 
and which generally characterise the simulation grid; the Independent parameters are those set by obser¬ 
vations; the Dependent column contains the parameters with an analytical dependency with respect to the 
dimensionless and/or independent parameters. 


Dimensionless 

Dependent 

Independent 

lO 

O 

T-1 

II 

£ 

(3 = avA 

B 0 = 10 14 G (observations) 

N z = 500 N y = 50 

™ 2 fiy/a 

va = 10 9 cm/s (correspond¬ 
ing to -B 0 ) 

AT = AZ = AY = 1 



xe (0,1,5,10,50,100) 



e = 0.3 



(7=1 
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where X s denote the numerical values of X measured at time s, and (} denotes averaging over all values s, 
respectively. Then the PSD is defined as [38] 

/ +oo 

Rx(r)e-^ fT dr. (27) 

-OO 

From a physical point of view the slope of the PSD of a time series X provides some statistical insight 
to the degree of correlation physical processes have with themselves. The observational light curves of the 
GRBs, as well as the simulated ones, are analysed in the next Section. 

3 Results and astrophysical implications 

We ran simulations in which we varied the time profile of the dimensionless velocity 14 according to the laws 

Vfe = constant, 14 ~ Vfc -1 ,14 ~ ft -1 , (28) 

and the value of y. These velocity profiles were chosen to mimic an ejected plasma flow placed in a gravi¬ 
tational field for a few specific cases: for 14 = constant the ejecta is so energetic that it does not feel any 
gravitational pull and the initial value of the velocity is conserved; the profiles 14 ~ Vk _1 and 14 ~ fc _1 are 
implemented to take into account a temporal decrease of velocity (equivalently, the plasma velocity decreases 
as a function of distance to the expulsion site). 

The main results of the simulations are the timestep evolution of the control parameter, the timestep 
evolution of the energy release, represented in the right panels of Figs. [T] - [3l the event lifetime distribution 
(FigQJ, the event size distribution (Fig. [5]) and the values of the spectral slope of the event lifetime distribution 
(Fig. [7]) and of the event size distribution (Fig. [8]), for different values of y. 

To allow a quick comparison with the observational data, in the left panels of Figs. [T]-El we present the 
light curves (counts) for three GRBs, with the plots done with the data file supplied by NASA’s HEASARC 
Data access database mmm- 

We identify the number of events needed to relax one critical onset with the lifetime of this avalanche. 
Thus, the distribution for lifetimes is equivalent to the distribution of number of events, 

D(N)~N~ aN . (29) 

For low advective velocities (y ~ 1) the exponent a at of a power-law distribution is approximately 1.6, which 
is in agreement with the theoretical value for two dimensional grids (see e.g. Eq. (3.7) from [3]). It is 
customary to argue that the exponent as of the power-law fitting the event size distribution is equal to the 
dimensionality of the system (i.e., as = 2 for a two-dimensional system) [39] . This is based on an unbiased 
diffusive random walk argument, where the distance L reached by a random walker in time T is L ~ \fT. 
More precisely, the walker is equally likely to choose any direction for his next step. This is not the case 
in our work, due to the inherent asymmetries of a deterministic background flow. As such, the numerical 
values of as depart from the value of 2. 

In order to produce a more quantitative comparison, we calculate the PSD of both the observed and 
simulated light curves; this is done by taking the corresponding time series, calculating their correlation 
functions, and then taking the Fourier transform. We then fit a model of the type PSD(f) ~ f~ a to them. 
The results from this fit, shown in Table [2] agree well with the results from the simulations as presented in 
Table [3] For illustration purposes, the PSD for GRB100919884 is shown on the same plot with the PSD of 
a simulated light curve in Fig. [6] The R 2 statistics presented in the third column of both these Tables is 
obtained following a linear regression scheme, based on minimizing the sum of the squared residuals; it is 
calculated as the ratio of the ratio of the model sum of squares to the total sum of squares. 

One important problem in CA-SOC simulations is the connection (conversion factor) between the timestep 
in the simulations and real time. This problem is even greater in simulations such as ours, which include 
two timescales: the advection timescale, set by the upwards propagating shock wave and the SOC timescale, 
i.e., the time needed by an avalanche to take place. 

As we see it, the main problem here is that, after SOC sets in, one timestep is the time needed for an 
avalanche to occur and this is not a constant. But in order to really be able to compare between observations 
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Figure 4: The lifetime distribution log(D) is plotted against the number of events log(IV), for x = 100 and 
different velocity profiles. 



log (S) 


Figure 5: The size distribution log (p(S)) is plotted against the size of the avalanche log(S'), for \ = 100 and 
different velocity profiles. 


Table 2: PSD of observations; a linear fit log PSD(f) = —alog/ was done. The third column shows the 
value of the R 2 statistics for the linear fit. _ 


Data identification 

a 

B' z 

GRB140919636 

0.708 

0.806 

GRB141223240 

0.585 

0.775 

GRB080805496 

0.746 

0.804 

GRB130702004 

0.960 

0.832 

GRB120119170 

1.505 

0.881 

GRB100919884 

0.912 

0.810 

GRB131230808 

0.722 

0.788 

GRB101126198 

1.698 

0.881 

GRB091209001 

1.006 

0.827 
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Table 3: PSD of simulations; a linear fit log PSD(f) = —a log/ was done. The third column shows the 
value of the R 2 statistics for the linear fit. _ 


Setup identification 

a 

H z 

- i /k, x = i 

1.125 

0.877 

v k - i/fc, x = so 

1.119 

0.890 

V k ~ l/Vk, x = i 

1.160 

0.880 

~ l/Vk, x = 10 

1.118 

0.876 

14 = const , x = 10 

1.703 

0.915 



log (f) 


Figure 6: Comparison of the variation of log(PSD) as a function of log(/) for the light curve of GRB 
GRB100919884 (upper curve) and for the simulated light curve with 14 ~ l/\/fc and x = 10 (lower curve). 


and simulations, a connection between real lapsed time and simulation timestep is proposed as follows: the 
magnetic reconnection rate in the Sun is of the order of ICC 2 — 10 - 3 eh si] ; reconnection rates in various 
astrophysical contexts do not vary that much, and hence we assume that this is also the order of magnitude 
of the reconnection in the GRB flow; we thus set a = 10 -4 ; we do this because a needs to be in a relationship 
to the smallest timescale in the problem such that all possible time-features can be probed by the simulation, 

i.e., OL T rec0 nnection- 

Next, we assume that once reconnection begins, the changing field line stresses will be transmitted 
throughout the simulation domain at the Alfven speed [25 j . Thus, the reconnection time is equal to the 
Alfven time, which is 10~ 4 s for our case (calculated using Table [T]). Then we assume that one timestep 
can be converted to seconds by multiplying the duration of a reconnection event by the average number of 
events produced by the simulation, estimated to be around 10 (see FigureS]). In the end, we conclude that 
10 3 timesteps make for one real second of observation time. 

In order to estimate the energy in dimensions of erg produced by the simulation we use Table [1] to get 
the order of magnitude for an E r produced by one timestep, E( lmestep ss 4 • 10 38 erg; if we assume that 
the conversion factor from time step to second is 10 3 , than the conversion factor between the dimensionless 
energy and the physical real energy is approximately 4 x 10 41 erg. The observed energy output in GRB 
explosions is around 10 44 erg [6]; the interval of dimensionless value for the energy obtained by simulations 
covers the values needed to produce agreement with observations. 

We note one very interesting fact: for zero velocity (i.e., no pulse at the base of the simulation grid), the 
model produces an event lifetime distribution slope in agreement with theory. For low impulse velocity, the 
event lifetime distribution slope is still close to the theoretic one. However, when x increases, for certain 
velocity profiles, the slope of D(N) vs. N departs from its theoretical value. The issue is best illustrated 
by plotting the value of the spectral slope as a function of x (Fig- 0- It is clear that for high values of 
X, one obtains a spectral slope with a value very close to 1, the theoretical accepted spectral slope for ID 
SOC. From a strict theoretical point of view, one dimensional magnetic reconnection is not possible; it was 
thus a problem (swept under the rug) to say that while the observational data for GRBs shows that a one 
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Table 4: The first column contains valid inquiries about observations, simulations and theory 
and how and if they agree. The second column contains an assessment of whether or not it 
is suitable to compare the simulation output with observational data. The third and fourth 
column contain the results of such comparison with observations and theory respectively. 
We use a — mark when such a comparison cannot be made, a •/ mark when the comparison 
produces good results and a x when the comparison fails. 


Output of the simulation 

May be compared with observational 
CRB data 

Obs. 

Theory 

Event lifetime distribution (an¬ 
alytical dependency N(E) ~ 
N~ a ) 

yes 

/ 

✓ 

PSD (analytical dependency 

p(f) ~ r a ) 

yes 

/ 

/ 

Light curve 

yes, a qualitative visual comparison 

/ 

- 

Maximum released energy 

yes, if the proposed conversion factor 
from dimensionless to dimensional is ac¬ 
cepted 

/ 


Number of peaks in the light 

curve 

yes, for some types of data 

/ 


Waiting time in the LC 

yes, if the waiting time is defined as 
time between subsequent avalanches in 
the grid (SOC characteristic) 
no, if waiting time is defined as lapsed 
time between initial pulse and maximum 
of the LC, because depends on the initial 
condition in the grid which is unknown 
no, if waiting time is defined as time 
lapsed between independent prominent 
light curve features, i.e. as those caused 
by different pulses at the base 


✓ 

Critical parameter timestep evo¬ 
lution 

no, this is just for SOC validation. For 
the moment this is just a theoretical tool. 


✓ 

Critical parameter maximum 
value 

no, although this could be done in prin¬ 
ciple if magnetograms of the studied re¬ 
gions would be available 



Magnetic field divergence 

no, this is just for general model valida¬ 
tion (fundamental physics constraint) 


✓ 


12 





































2.5 


2.0 


1.5 


1.0 


10 


15 20 

*/io 


25 30 


Figure 7: A fit of the type D(N) ~ A" is made for the lifetime distribution. The plot shows the values 
of oat as a function of %, for V ~ Vfe -1 . 
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Figure 8: A fit of the type p(S) ~ 5 “ s is made for the event size distribution. The plot shows the values 
of as as a function of %, for R ~ \/A; -1 . 


dimensional SOC is occurring , this is due to magnetic reconnection. 

However, our results can be interpreted as showing that the physics is genuinely 2D, but an energetic 
initial impulse leads to a ID SOC signature. 

The physics and observations of GRBs are extremely complex. It is understandable that a simple CA 
model cannot and does not reproduce all features of the real process. As such, we list some of the caveats 
of our model: it does not produce information about bulk Lorentz factor, about isotropy (or lack thereof) 
of the radiation. Also, although under some assumptions, magnetic CA models can produce maps of the 
magnetic field on the microscopic level, this was not the purpose of this paper. A discussion of how well the 
model agrees with observational data is presented in Table 0] 


4 Conclusions 

In the present paper we have presented a two dimensional CA model for a magnetized grid with an initial large 
flow and a stochastic perturbation. The discretized laws for the evolution of the magnetic field were deduced 
and implemented. This study extends to two dimensions the one dimensional SOC approach introduced 
in |18] for the analysis of one dimensional magnetized plasma flows. As a main result of our investigation 
we have found that the two dimensional magnetized plasma with a background flow system exhibits SOC. 
The event lifetime and event size distributions and the emitted electromagnetic energy were obtained, and 
their astrophysical interpretation was briefly discussed. It was shown that this approach can be used to 
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model GRBs and some of their properties, as summarized in Tabic [I] Most importantly, the model produces 
results which agree with the observations on the same time and energy scales, which to our knowledge is 
a novelty. [55] suggested to use CA and SOC methods to study and interpret the X-ray flares of GRBs. 
Here we further postulated that the main burst, i.e., the GRB itself, also arise from SOC processes. In 
contrast to the standard GRB models, where each pulse correspond to one energy injection (e.g., internal 
shock model), here we show that in the SOC model one energy input is able to generate a series of pulses, 
which are characteristic of many GRBs. Moreover, we have found that some general features of GRBs (and 
especially the light curves) can be described in the framework of the present SOC model. 

It would be interesting to compare our mathematical SOC formalism with the one used to study the 
properties of the accretion by black holes. In the model introduced in [55115U] . the accretion disk is divided 
into two parts, the outer disk region, and the inner disk region, respectively. In the outer disk region, the 
gas drifts smoothly inward, while in the inner disk region blobs are formed. The whole disk is divided into 
64 rings, and angularly divided into 64 equal parts. The physical quantity stored in the cells is the mass in 
that region. Since the mass is continuously flowing to the outermost region of the disk, a randomly chosen 
cell in the ring receives a fixed amount of mass m. If the mass of the cell is larger than the critical value 
M cr it, a mass flow will occur to the three nearest cells in the adjacent inner ring. By denoting the mass in 
the cell with coordinates i (representing the radial position) and j (representing the angular position) by 
Mjj, when Mij > M cr n, the following CA process happens: Mjj —> Mij — 3m, M ,;_ij —> Mj_ij + m, 
Mi- i,j±i —> Mj_i j±i + m. The total X-ray luminosity is of the order of the gravitational potential energy 
lost in the process, and it can be obtained as L ss ^ (GMm/ri-i — GMm/ri). 

In contrast to the above approach, in the present paper we start from the basic result that diffusion type 
partial differential equations driven by noise can produce characteristic effects similar to SOC [35]. Hence, 
we consider a discretized version of the magnetic induction equation, in the presence of a stochastic term 
generated by the matter flow and the gradient of the magnetic field. The magnetic induction equation is 
considered in two limiting cases, corresponding to the diffusive and advective regimes, respectively. There are 
also significant differences in the physical approach of the two models. While in the accretion disk model of 
mm the parameter stored in the grid is the mass, in the present approach the stored quantity is the value 
of the magnetic field. The disk luminosity is assumed to be purely of gravitational origin in the accretion 
disk CA models, while in the present case the radiation emission is purely electromagnetic. 

An important characteristic of systems that reached SOC is that a slight change in the external conditions 
of the system determines the generation of avalanches of various sizes. These avalanches keep, on average, 
the system near the critical state [3 Oil- It is important to note that a minor effect can result in an avalanche 
of a much bigger size (catastrophe). This kind of catastrophic behaviour could be also responsible for the 
GRBs explosions, and overall evolution. Therefore SOC and CA methods can be used, and may prove useful, 
in modelling early phases - explosive or immediately post-explosive - of the GRBs. Once a given critical 
condition is reached, it is self-maintained, i.e., there is no need to vary or fine tune its physical and control 
parameters. This situation is quite different as compared, for example, with the critical state in liquid, where 
two thermodynamic parameters(the temperature and the density) must be tuned to keep the fluid system 
near its critical point. 

It is a general property of SOC that the event lifetime and size distributions of avalanches are described 
in terms of power laws with negative exponents. This is one of the main reasons why models based on SOC 
can be successfully used to describe various natural phenomena. In the present paper we have developed 
some necessary (but not sufficient) tools that are required for the in depth comparison of the theoretical 
models with the specific observations of the cosmological Gamma Ray Bursts. 
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